Exploring covariates, rs-fMRI usability, and functional connectivity
Liwei Wang and MB Nebel
9/10/2021
#for flowchart
library(DiagrammeR)
#to export flowchart as .png
#webshot::install_phantomjs()
library(DiagrammeRsvg)
library(rsvg)
library(grid)
library(xtable)
library(dplyr)
library(ggplot2)
library(gridExtra)
#nest/unnest
library(tidyr)
#map function (kind of like a for loop)
library(purrr)
#tidy model summary
library(broom)
library(readxl)
#for tableby
library(arsenal)
#Also uses functions from plyr, scales, mgcv, cowplotDefine geom_split_violin()
Based on https://debruine.github.io/posts/plot-comparison/
GeomSplitViolin <- ggproto("GeomSplitViolin", GeomViolin, draw_group = function(self,
data, ..., draw_quantiles = NULL) {
data <- transform(data, xminv = x - violinwidth * (x - xmin), xmaxv = x + violinwidth *
(xmax - x))
grp <- data[1, "group"]
newdata <- plyr::arrange(transform(data, x = if (grp%%2 == 1)
xminv else xmaxv), if (grp%%2 == 1)
y else -y)
newdata <- rbind(newdata[1, ], newdata, newdata[nrow(newdata), ], newdata[1,
])
newdata[c(1, nrow(newdata) - 1, nrow(newdata)), "x"] <- round(newdata[1, "x"])
if (length(draw_quantiles) > 0 & !scales::zero_range(range(data$y))) {
stopifnot(all(draw_quantiles >= 0), all(draw_quantiles <= 1))
quantiles <- ggplot2:::create_quantile_segment_frame(data, draw_quantiles)
aesthetics <- data[rep(1, nrow(quantiles)), setdiff(names(data), c("x", "y")),
drop = FALSE]
aesthetics$alpha <- rep(1, nrow(quantiles))
both <- cbind(quantiles, aesthetics)
quantile_grob <- GeomPath$draw_panel(both, ...)
ggplot2:::ggname("geom_split_violin", grid::grobTree(GeomPolygon$draw_panel(newdata,
...), quantile_grob))
} else {
ggplot2:::ggname("geom_split_violin", GeomPolygon$draw_panel(newdata, ...))
}
})
geom_split_violin <- function(mapping = NULL, data = NULL, stat = "ydensity", position = "identity",
..., draw_quantiles = NULL, trim = TRUE, scale = "area", na.rm = FALSE, show.legend = NA,
inherit.aes = TRUE) {
layer(data = data, mapping = mapping, stat = stat, geom = GeomSplitViolin, position = position,
show.legend = show.legend, inherit.aes = inherit.aes, params = list(trim = trim,
scale = scale, draw_quantiles = draw_quantiles, na.rm = na.rm, ...))
}Load Data
load('../Data/DataWithPropensities_seed1.RData')
#convert PrimaryDiagnosis to factor
dat3$PrimaryDiagnosis <- factor(dat3$PrimaryDiagnosis, levels = c("Autism", "None"))
levels(dat3$PrimaryDiagnosis)## [1] "Autism" "None"dat3$PrimaryDiagnosis <- relevel(dat3$PrimaryDiagnosis, "None")
levels(dat3$PrimaryDiagnosis) = c("TD", "ASD")49 ASD-only + 353 TD = 402 ADHD_Secondary=0
Reshape data to combine motion QC levels
# create dummy factor to include all subjects
dat3$noExclusion <- ifelse(dat3$ID > 0, "Pass", "Pass")
dat3$noExclusion <- factor(dat3$noExclusion, levels = c("Pass", "Fail"))
# convert KKI_criteria to factor with reference level 'Pass'
dat3$KKI_criteria <- factor(dat3$KKI_criteria, levels = c("Pass", "Fail"))
# convert Ciric_length to factor with reference level 'Pass' to match
# KKI_criteria
dat3$Ciric_length <- factor(dat3$Ciric_length, levels = c("Pass", "Fail"))
# combine Ciric_length, KKI, and noExclusion exclusion into one variable
allVariables = c("ID", "PrimaryDiagnosis", "AgeAtScan", "Ciric_length", "KKI_criteria",
"noExclusion", "PANESS.TotalOverflowNotAccountingForAge", "SRS.Score", "WISC.GAI",
"DuPaulHome.InattentionRaw", "DuPaulHome.HyperactivityRaw", "ADOS.Comparable.Total",
"CurrentlyOnStimulants", "HeadCoil", "Sex", "ADHD_Secondary", "SES.Family", "Race2",
"handedness", "CompletePredictorCases", "YearOfScan", "MeanFramewiseDisplacement.KKI")
idVariables = c("ID", "PrimaryDiagnosis", "AgeAtScan", "PANESS.TotalOverflowNotAccountingForAge",
"SRS.Score", "WISC.GAI", "DuPaulHome.InattentionRaw", "DuPaulHome.HyperactivityRaw",
"ADOS.Comparable.Total", "CurrentlyOnStimulants", "HeadCoil", "Sex", "ADHD_Secondary",
"SES.Family", "Race2", "handedness", "CompletePredictorCases", "YearOfScan",
"MeanFramewiseDisplacement.KKI")
qcMelt <- reshape2::melt(dat3[, allVariables], id.vars = names(dat3)[which(names(dat3) %in%
idVariables)], variable.name = "Motion.Exclusion.Level", value.name = "Included")
# rename exclusion levels NOTE: need None to be highest level for
# geom_split_violin
levels(qcMelt$Motion.Exclusion.Level) <- c("Strict", "Lenient", "None")
# convert Included to factor with pass as reference
qcMelt$Included <- factor(qcMelt$Included, levels = c("Pass", "Fail"))
# rename levels of value
levels(qcMelt$Included) <- c("Included", "Excluded")
with(qcMelt, table(PrimaryDiagnosis, Motion.Exclusion.Level, Included))## , , Included = Included
##
## Motion.Exclusion.Level
## PrimaryDiagnosis Strict Lenient None
## TD 151 308 372
## ASD 29 114 173
##
## , , Included = Excluded
##
## Motion.Exclusion.Level
## PrimaryDiagnosis Strict Lenient None
## TD 221 64 0
## ASD 144 59 0Lenient motion QC = KKI_criteria
Strict motion QC = Ciric_length
Limit initial data set to complete cases
dat3 <- filter(dat3, CompletePredictorCases==1)1. Impact of motion QC on sample size
completeCases <- filter(qcMelt, CompletePredictorCases==1)
tabN<- tableby(Motion.Exclusion.Level~
Included, data=filter(completeCases, Motion.Exclusion.Level!="None"))
summary(tabN,
title='Proportion of complete cases included/excluded by motion QC',digits=0, digits.p=4,digits.pct=1, numeric.simplify=TRUE, total=FALSE, test=FALSE)| Strict (N=504) | Lenient (N=504) | |
|---|---|---|
| Included | ||
| Included | 171 (33.9%) | 403 (80.0%) |
| Excluded | 333 (66.1%) | 101 (20.0%) |
tabASD<- tableby(Motion.Exclusion.Level~Included, data=filter(completeCases,
PrimaryDiagnosis=="ASD"))
summary(tabASD,
title = "Proportion of ASD complete cases included/excluded",
digits=0, digits.p=4,digits.pct=1, numeric.simplify=TRUE, total=FALSE, test=FALSE)| Strict (N=151) | Lenient (N=151) | None (N=151) | |
|---|---|---|---|
| Included | |||
| Included | 29 (19.2%) | 107 (70.9%) | 151 (100.0%) |
| Excluded | 122 (80.8%) | 44 (29.1%) | 0 (0.0%) |
tabTD<- tableby(Motion.Exclusion.Level~Included, data=filter(completeCases,
PrimaryDiagnosis=="TD"))
summary(tabTD,
title = "Proportion of TD complete cases included/excluded",
digits=0, digits.p=4,digits.pct=1, numeric.simplify=TRUE, total=FALSE, test=FALSE)| Strict (N=353) | Lenient (N=353) | None (N=353) | |
|---|---|---|---|
| Included | |||
| Included | 142 (40.2%) | 296 (83.9%) | 353 (100.0%) |
| Excluded | 211 (59.8%) | 57 (16.1%) | 0 (0.0%) |
1b. Summarize sociodemographic info for complete cases for paper
#make M reference level for sex
completeCases$Sex <- relevel(as.factor(completeCases$Sex), "M")
#use chisq for Sex and handedness, kruskal-wallis rank test for Age
tabSex<- tableby( PrimaryDiagnosis~Sex+AgeAtScan+handedness, data=filter(completeCases, Motion.Exclusion.Level=="None" & Included=="Included"), control=tableby.control(numeric.test="kwt", cat.test="chisq", total=FALSE))
#use fisher's exact test for Race, k-w for SES
tabRace<- tableby( PrimaryDiagnosis~Race2+SES.Family, data=filter(completeCases, Motion.Exclusion.Level=="None" & Included=="Included"), control=tableby.control(numeric.test="kwt", cat.test="fe", total=FALSE))
tab12 <- merge(tabSex, tabRace)
#Currently on Stimulants - no test because no TDs are currently on stimulants by design
completeCases$CurrentlyOnStimulants <- factor(completeCases$CurrentlyOnStimulants)
#rename factor levels
levels(completeCases$CurrentlyOnStimulants)[levels(completeCases$CurrentlyOnStimulants)=="0"] <- "No"
levels(completeCases$CurrentlyOnStimulants)[levels(completeCases$CurrentlyOnStimulants)=="1"] <- "Yes"
tabStim<- tableby( PrimaryDiagnosis~CurrentlyOnStimulants, data=filter(completeCases, Motion.Exclusion.Level=="None" & Included=="Included"), control=tableby.control(total=FALSE, test=FALSE))
tab123 <- merge(tab12, tabStim)
summary(tab123)| TD (N=353) | ASD (N=151) | p value | |
|---|---|---|---|
| Sex | < 0.001 | ||
| M | 245 (69.4%) | 127 (84.1%) | |
| F | 108 (30.6%) | 24 (15.9%) | |
| AgeAtScan | 0.826 | ||
| Mean (SD) | 10.363 (1.248) | 10.324 (1.363) | |
| Range | 8.020 - 12.980 | 8.010 - 12.990 | |
| handedness | 0.253 | ||
| Left | 17 (4.8%) | 12 (7.9%) | |
| Mixed | 19 (5.4%) | 11 (7.3%) | |
| Right | 317 (89.8%) | 128 (84.8%) | |
| Race2 | 0.004 | ||
| African American | 36 (10.2%) | 9 (6.0%) | |
| Asian | 27 (7.6%) | 3 (2.0%) | |
| Biracial | 45 (12.7%) | 12 (7.9%) | |
| Caucasian | 245 (69.4%) | 127 (84.1%) | |
| SES.Family | 0.006 | ||
| Mean (SD) | 54.135 (9.390) | 51.964 (9.379) | |
| Range | 18.500 - 66.000 | 27.000 - 66.000 | |
| CurrentlyOnStimulants | |||
| No | 353 (100.0%) | 97 (64.2%) | |
| Yes | 0 (0.0%) | 54 (35.8%) |
Generate version of table to paste into overleaf
tab <- xtable(as.data.frame(summary(tab123)))
print(tab, type="latex")## % latex table generated in R 4.1.1 by xtable 1.8-4 package
## % Thu Oct 7 00:43:14 2021
## \begin{table}[ht]
## \centering
## \begin{tabular}{rllll}
## \hline
## & & TD (N=353) & ASD (N=151) & p value \\
## \hline
## 1 & **Sex** & & & $<$ 0.001 \\
## 2 & \ \ \ M & 245 (69.4\%) & 127 (84.1\%) & \\
## 3 & \ \ \ F & 108 (30.6\%) & 24 (15.9\%) & \\
## 4 & **AgeAtScan** & & & 0.826 \\
## 5 & \ \ \ Mean (SD) & 10.363 (1.248) & 10.324 (1.363) & \\
## 6 & \ \ \ Range & 8.020 - 12.980 & 8.010 - 12.990 & \\
## 7 & **handedness** & & & 0.253 \\
## 8 & \ \ \ Left & 17 (4.8\%) & 12 (7.9\%) & \\
## 9 & \ \ \ Mixed & 19 (5.4\%) & 11 (7.3\%) & \\
## 10 & \ \ \ Right & 317 (89.8\%) & 128 (84.8\%) & \\
## 11 & **Race2** & & & 0.004 \\
## 12 & \ \ \ African American & 36 (10.2\%) & 9 (6.0\%) & \\
## 13 & \ \ \ Asian & 27 (7.6\%) & 3 (2.0\%) & \\
## 14 & \ \ \ Biracial & 45 (12.7\%) & 12 (7.9\%) & \\
## 15 & \ \ \ Caucasian & 245 (69.4\%) & 127 (84.1\%) & \\
## 16 & **SES.Family** & & & 0.006 \\
## 17 & \ \ \ Mean (SD) & 54.135 (9.390) & 51.964 (9.379) & \\
## 18 & \ \ \ Range & 18.500 - 66.000 & 27.000 - 66.000 & \\
## 19 & **CurrentlyOnStimulants** & & & \\
## 20 & \ \ \ No & 353 (100.0\%) & 97 (64.2\%) & \\
## 21 & \ \ \ Yes & 0 (0.0\%) & 54 (35.8\%) & \\
## \hline
## \end{tabular}
## \end{table}1c. Create exclusion flowchart
1d. Proportion excluded
Define theme for proportion excluded plots
#Easier to make three separate panels and combine them than to use faceting function
My_Theme_prop = theme_light()+theme(
legend.title =element_blank(),
axis.title.x = element_text(size = 12),
axis.title.y = element_text(size = 10),
plot.title = element_text(size = 30),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 10),
strip.text.x = element_text(size = 10,color="black"),
strip.background = element_rect(fill="white"))1d.1 Proportion excluded in each Diagnostic Group
motion <- filter(completeCases, Motion.Exclusion.Level != "None")
motion$Motion.Exclusion.Level <- droplevels(motion$Motion.Exclusion.Level)
motion <- group_by(motion, PrimaryDiagnosis, Motion.Exclusion.Level, Included)
dx_proportions <- ggplot(motion, aes(x = PrimaryDiagnosis, fill = Included)) + geom_bar(position = "fill",
alpha = 0.6) + facet_grid(~Motion.Exclusion.Level) + scale_fill_manual(values = c("#FDE599",
"#9FB0CC")) + scale_color_manual(values = c("#E9D38D", "#8C9AB4")) + My_Theme_prop +
theme(legend.title = element_blank()) + theme(legend.title = element_blank()) +
ylab("Proportion of Children") + theme(legend.position = "right") + theme(legend.margin = margin(t = 0,
r = 0, b = -1, l = -1)) + theme(legend.key.size = unit(0.08, "in")) + theme(axis.title.x = element_blank()) +
theme(axis.text.x = element_text(size = 10))
png("Figures/fig_propExcludedDx_cc.png", width = 3, height = 2, units = "in", res = 200)
dx_proportions
invisible(dev.off())
dx_proportions
dx_proportions <- dx_proportions + theme(legend.position = "none") + theme(axis.text.x = element_text(size = 8),
axis.title.x = element_text(size = 12)) + xlab("Diagnosis Group")
# Pearson's chi squared tests
extib <- tibble(motion)
exNest <- extib %>%
select(c("PrimaryDiagnosis", "Motion.Exclusion.Level", "Included")) %>%
group_by(Motion.Exclusion.Level) %>%
tidyr::nest()
# nested models
ex_chisq <- exNest %>%
mutate(stats = map(data, ~broom::tidy(chisq.test(.x$PrimaryDiagnosis, .x$Included)))) %>%
unnest(stats)
ex_chisq## # A tibble: 2 × 6
## # Groups: Motion.Exclusion.Level [2]
## Motion.Exclusion.Level data statistic p.value parameter method
## <fct> <list> <dbl> <dbl> <int> <chr>
## 1 Strict <tibble [504 × 2]> 19.9 8.07e-6 1 Pearso…
## 2 Lenient <tibble [504 × 2]> 10.3 1.30e-3 1 Pearso…The proportion of children excluded differs across diagnostic groups using both the lenient and strict motion QC.
1d.2 Currently on Stimulants
stim <- filter(motion, PrimaryDiagnosis=="ASD")
stim$PrimaryDiagnosis <- droplevels(stim$PrimaryDiagnosis)
stim$CurrentlyOnStimulants <- factor(stim$CurrentlyOnStimulants)
levels(stim$CurrentlyOnStimulants)## [1] "No" "Yes"levels(stim$CurrentlyOnStimulants) = c("No", "Yes")
stimNest <- tibble(stim) %>%
select(c("CurrentlyOnStimulants",
"Motion.Exclusion.Level", "Included")) %>%
group_by(Motion.Exclusion.Level) %>%
tidyr::nest()
#nested models
stim_chisq <- stimNest %>%
mutate(stats = map(data, ~broom::tidy(chisq.test(.x$CurrentlyOnStimulants, .x$Included)))) %>%
unnest(stats)
stim_chisq## # A tibble: 2 × 6
## # Groups: Motion.Exclusion.Level [2]
## Motion.Exclusion.Level data statistic p.value parameter method
## <fct> <list> <dbl> <dbl> <int> <chr>
## 1 Strict <tibble [151 × 2]> 6.40 0.0114 1 Pearson…
## 2 Lenient <tibble [151 × 2]> 6.39 0.0115 1 Pearson…1d.5 Proportion excluded by Currently on stimulants in ASD group
stim_proportions <- ggplot(stim, aes(x = CurrentlyOnStimulants, fill = Included))+
geom_bar(position = "fill", alpha = .6)+
facet_grid(~Motion.Exclusion.Level)+
scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
My_Theme_prop+
xlab("Currently On Stimulants?")+
theme(legend.title = element_blank())+
#theme(legend.position="top")+
theme(legend.text = element_text(size=7))+
ylab("Proportion of Autistic Children")+
theme(legend.position = "none")+
theme(strip.text.x = element_blank())
#theme(plot.margin = unit(c(1, 6, 1, 6), "cm"))
prop_legend = cowplot::get_legend(stim_proportions + guides(color = guide_legend(ncol = 1))+
theme(legend.position = "bottom",
legend.text = element_text(size = 11),
legend.key.size=unit(.1, "in")))1d.6 Combine proportion plot with legend & save
####Impact of motion QC on framewise displacement metrics
#mean and max FD are different for lenient and strict because some frames at the beginning and/or end of the scan are excluded for scans to pass lenient motion QC
dat3$MeanFD.None = dat3$MeanFramewiseDisplacement
dat3$MaxFD.None = dat3$MaxFramewiseDisplacement
#same for all levels
dat3$FramesWithFDLessThanOrEqualTo250microns.None = dat3$FramesWithFDLessThanOrEqualTo250microns
dat3$FramesWithFDLessThanOrEqualTo250microns.KKI = dat3$FramesWithFDLessThanOrEqualTo250microns
meanFD = c("ID", "MeanFramewiseDisplacement", "MeanFramewiseDisplacement.KKI",
"MeanFD.None")
maxFD = c("ID", "MaxFramewiseDisplacement", "MaxFramewiseDisplacement.KKI", "MaxFD.None")
framesFD = c("ID", "FramesWithFDLessThanOrEqualTo250microns",
"FramesWithFDLessThanOrEqualTo250microns.KKI",
"FramesWithFDLessThanOrEqualTo250microns.None")
fdID = c("ID")
meanFD.df <- reshape2::melt(dat3[, meanFD],
id.vars=names(dat3)[which(names(dat3) %in% fdID)],
variable.name = "Motion.Exclusion.Level",
value.name = "MeanFramewiseDisplacement")
levels(meanFD.df$Motion.Exclusion.Level)## [1] "MeanFramewiseDisplacement" "MeanFramewiseDisplacement.KKI"
## [3] "MeanFD.None"#rename levels to match motion QC levels in completeCases
levels(meanFD.df$Motion.Exclusion.Level) <- c("Strict", "Lenient", "None")
#repeat for MaxFD
maxFD.df <- reshape2::melt(dat3[, maxFD],
id.vars=names(dat3)[which(names(dat3) %in% fdID)],
variable.name = "Motion.Exclusion.Level",
value.name = "MaxFramewiseDisplacement")
levels(maxFD.df$Motion.Exclusion.Level)## [1] "MaxFramewiseDisplacement" "MaxFramewiseDisplacement.KKI"
## [3] "MaxFD.None"#rename levels to match motion QC levels in completeCases
levels(maxFD.df$Motion.Exclusion.Level) <- c("Strict", "Lenient", "None")
#merge meanFD.df and maxFD.df
fdMerg <- merge(meanFD.df, maxFD.df)
#repeat for FramesWithFDLessThanOrEqualTo250microns
frames.df <- reshape2::melt(dat3[, framesFD],
id.vars=names(dat3)[which(names(dat3) %in% fdID)],
variable.name = "Motion.Exclusion.Level",
value.name = "FramesWithFDLessThanOrEqualTo250microns")
levels(frames.df$Motion.Exclusion.Level)## [1] "FramesWithFDLessThanOrEqualTo250microns"
## [2] "FramesWithFDLessThanOrEqualTo250microns.KKI"
## [3] "FramesWithFDLessThanOrEqualTo250microns.None"#rename levels to match motion QC levels in qcMelt
levels(frames.df$Motion.Exclusion.Level) <- c("Strict", "Lenient", "None")
#merge with fdMerg
fdMerg <- merge(fdMerg, frames.df)
#merge with completeCases
completeCases <- merge(completeCases, fdMerg)
passOnly <- filter(completeCases, Included=="Included")mean framewise displacement split violin none, lenient, strict
meanFD_violin <- ggplot(passOnly, aes(Motion.Exclusion.Level, MeanFramewiseDisplacement,
fill = PrimaryDiagnosis, color=PrimaryDiagnosis)) +
geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE, adjust = 2) +
stat_summary(fun = "mean", position = position_dodge(width = 0.18),
color="black", geom="point", aes(y=MeanFramewiseDisplacement))+
stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.18),
color="black", geom="errorbar", width=.15)+
scale_fill_manual(values = c("#009E73", "#FDE599"))+
scale_color_manual(values = c("#05634a", "#E9D38D"))+
My_Theme_prop+theme(legend.title =element_blank())+
theme(axis.text.y=element_text(size=8))+
theme(legend.text = element_text(size = 7))+
theme(axis.title.y = element_blank())+
theme(axis.title.x = element_blank())+
theme(legend.position = c(0.35, .7))+
theme(legend.key.size=unit(.15, "in"))+
theme(legend.box.margin = margin(-2, -2, -2, -2))+
ggtitle("Mean FD")+
theme(plot.title = element_text(size = 8, hjust=0.5))max framewise displacement split violin none, lenient, strict
maxFD_violin <- ggplot(passOnly, aes(Motion.Exclusion.Level, MaxFramewiseDisplacement,
fill = PrimaryDiagnosis, color=PrimaryDiagnosis)) +
geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE, adjust = 2) +
stat_summary(fun = "mean", position = position_dodge(width = 0.18), color="black",
geom="point", aes(y=MaxFramewiseDisplacement))+
stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.18),
color="black", geom="errorbar", width=.15)+
scale_fill_manual(values = c("#009E73", "#FDE599"))+
scale_color_manual(values = c("#05634a", "#E9D38D"))+
My_Theme_prop+
theme(axis.text.y=element_text(size=8))+
theme(axis.title.y = element_blank())+
theme(axis.title.x = element_blank())+
theme(legend.position = "none")+
ggtitle("Max FD")+
theme(plot.title = element_text(size = 8, hjust=0.5))frames with FD <= 0.25 mm split violin none, lenient, strict
frames_violin <- ggplot(passOnly, aes(Motion.Exclusion.Level, FramesWithFDLessThanOrEqualTo250microns,
fill = PrimaryDiagnosis, color=PrimaryDiagnosis)) +
geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE, adjust = 2) +
stat_summary(fun = "mean", position = position_dodge(width = 0.18), color="black",
geom="point", aes(y=FramesWithFDLessThanOrEqualTo250microns))+
stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.18),
color="black", geom="errorbar", width=.15)+
scale_fill_manual(values = c("#009E73", "#FDE599"))+
scale_color_manual(values = c("#05634a", "#E9D38D"))+
My_Theme_prop+
theme(axis.text.y=element_text(size=8))+
theme(axis.title.y = element_blank())+
theme(axis.title.x = element_blank())+
theme(legend.position = "none")+
ggtitle("Frames with FD<0.25 mm")+
theme(plot.title = element_text(size = 8, hjust=0.5))Combine 3 FD plots into one figure and save
## quartz_off_screen
## 2
Framewise displacement metrics are more similar across diagnostic groups using the strict motion QC, but very few participants are labeled as usable.
Filter out “None” from Motion.Exclusion.Level to make remaining group differences in FD metrics following motion QC easier to see
mean framewise displacement
passOnly <- filter(passOnly, Motion.Exclusion.Level!="None")
meanFD_violin <- ggplot(passOnly, aes(Motion.Exclusion.Level, MeanFramewiseDisplacement,
fill = PrimaryDiagnosis, color=PrimaryDiagnosis)) +
geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE, adjust = 2) +
stat_summary(fun = "mean", position = position_dodge(width = 0.18),
color="black", geom="point", aes(y=MeanFramewiseDisplacement))+
stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.18),
color="black", geom="errorbar", width=.15)+
scale_fill_manual(values = c("#009E73", "#FDE599"))+
scale_color_manual(values = c("#05634a", "#E9D38D"))+
My_Theme_prop+theme(legend.title =element_blank())+
theme(axis.text.y=element_text(size=8))+
theme(legend.text = element_text(size = 7))+
theme(axis.title.y = element_blank())+
theme(axis.title.x = element_blank())+
theme(legend.position = c(0.35, .7))+
theme(legend.key.size=unit(.15, "in"))+
theme(legend.box.margin = margin(-2, -2, -2, -2))+
ggtitle("Mean FD")+
theme(plot.title = element_text(size = 8, hjust=0.5))max framewise displacement
maxFD_violin <- ggplot(passOnly, aes(Motion.Exclusion.Level, MaxFramewiseDisplacement,
fill = PrimaryDiagnosis, color=PrimaryDiagnosis)) +
geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE, adjust = 2) +
stat_summary(fun = "mean", position = position_dodge(width = 0.18), color="black",
geom="point", aes(y=MaxFramewiseDisplacement))+
stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.18),
color="black", geom="errorbar", width=.15)+
scale_fill_manual(values = c("#009E73", "#FDE599"))+
scale_color_manual(values = c("#05634a", "#E9D38D"))+
My_Theme_prop+
theme(axis.text.y=element_text(size=8))+
theme(axis.title.y = element_blank())+
theme(axis.title.x = element_blank())+
theme(legend.position = "none")+
ggtitle("Max FD")+
theme(plot.title = element_text(size = 8, hjust=0.5))frames with FD <= 0.25 mm
frames_violin <- ggplot(passOnly, aes(Motion.Exclusion.Level, FramesWithFDLessThanOrEqualTo250microns,
fill = PrimaryDiagnosis, color=PrimaryDiagnosis)) +
geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE, adjust = 2) +
stat_summary(fun = "mean", position = position_dodge(width = 0.18), color="black",
geom="point", aes(y=FramesWithFDLessThanOrEqualTo250microns))+
stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.18),
color="black", geom="errorbar", width=.15)+
scale_fill_manual(values = c("#009E73", "#FDE599"))+
scale_color_manual(values = c("#05634a", "#E9D38D"))+
My_Theme_prop+
theme(axis.text.y=element_text(size=8))+
theme(axis.title.y = element_blank())+
theme(axis.title.x = element_blank())+
theme(legend.position = "none")+
ggtitle("Frames with FD<0.25 mm")+
theme(plot.title = element_text(size = 8, hjust=0.5))Combine 3 FD plots with just lenient and strict Motion.Exclusion.Levels and save
## quartz_off_screen
## 2
2. rs-fMRI exclusion probability as a function of age and symptom severity
Covariates
phenoVariables <- c("ID", "PrimaryDiagnosis",
"ADOS.Comparable.Total",
"SRS.Score",
"PANESS.TotalOverflowNotAccountingForAge",
"DuPaulHome.InattentionRaw",
"DuPaulHome.HyperactivityRaw",
"AgeAtScan",
"WISC.GAI",
"Motion.Exclusion.Level", "Included")
phenoIDs <- c("ID", "PrimaryDiagnosis", "Motion.Exclusion.Level", "Included")
aim1 <- reshape2::melt(completeCases[, phenoVariables],
id.vars=names(completeCases)[which(names(completeCases) %in% phenoIDs)])
levels(aim1$variable) <- c("ADOS", "SRS", "Motor Overflow", "Inattention",
"Hyperactivity", "Age", "GAI")
aim1G <- group_by(aim1, PrimaryDiagnosis, Motion.Exclusion.Level, Included, variable)2a. Univarate GAMs
aim1$delta = rep(NA,length=nrow(aim1))
aim1$delta = ifelse(aim1$Included=="Included",1,0)
aim1tib <- tibble(filter(aim1, Motion.Exclusion.Level!="None"))
aim1tib$Motion.Exclusion.Level <- droplevels(aim1tib$Motion.Exclusion.Level)
aim1Nest <- aim1tib %>%
group_by(Motion.Exclusion.Level, variable) %>%
tidyr::nest()
#nested models
nested_gams <- aim1Nest %>%
mutate(model = map(data, ~mgcv::gam(1-delta~s(value, k=-1), data = na.omit(.x),
family=binomial(link=logit), method="REML")),
coefs = map(model, tidy, conf.int = FALSE),
Rsq = map_dbl(model, ~summary(.)$r.sq)) %>%
unnest(coefs)
#Ben: correct for 7 lenient and 7 strict
nested_gams_len <- nested_gams %>%
filter(Motion.Exclusion.Level=="Lenient")
nested_gams_len$p.fdr = p.adjust(nested_gams_len$p.value, method = "BH")
nested_gams_strict <- nested_gams %>%
filter(Motion.Exclusion.Level=="Strict")
nested_gams_strict$p.fdr = p.adjust(nested_gams_strict$p.value, method = "BH")
#combine to print
nested_gams <- rbind(nested_gams_len, nested_gams_strict)
#list adjusted p values
nested_gams[, c(1:2,6:11)]## # A tibble: 14 × 8
## # Groups: Motion.Exclusion.Level, variable [14]
## Motion.Exclusion.… variable edf ref.df statistic p.value Rsq p.fdr
## <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Lenient ADOS 2.04 2.47 18.4 2.08e-4 0.0442 4.86e-4
## 2 Lenient SRS 1.83 2.29 17.4 2.95e-4 0.0466 5.17e-4
## 3 Lenient Motor Ove… 1.00 1.00 15.6 7.58e-5 0.0308 4.54e-4
## 4 Lenient Inattenti… 1.38 1.67 7.44 1.17e-2 0.0157 1.17e-2
## 5 Lenient Hyperacti… 2.15 2.70 13.2 4.59e-3 0.0242 5.36e-3
## 6 Lenient Age 1.00 1.00 9.33 2.26e-3 0.0164 3.17e-3
## 7 Lenient GAI 1.00 1.00 14.7 1.30e-4 0.0288 4.54e-4
## 8 Strict ADOS 1.00 1.00 21.9 2.79e-6 0.0444 9.76e-6
## 9 Strict SRS 1.00 1.00 22.7 1.53e-6 0.0567 9.76e-6
## 10 Strict Motor Ove… 1.00 1.00 10.2 1.42e-3 0.0194 1.99e-3
## 11 Strict Inattenti… 1.00 1.00 17.8 2.47e-5 0.0360 4.31e-5
## 12 Strict Hyperacti… 1.88 2.35 25.0 1.30e-5 0.0516 3.04e-5
## 13 Strict Age 1.76 2.20 7.73 2.84e-2 0.0129 2.84e-2
## 14 Strict GAI 1.00 1.00 7.55 6.02e-3 0.0130 7.02e-3#max p value for 7 lenient models
max(nested_gams_len$p.fdr)## [1] 0.01171046#max p value for 7 strict models
max(nested_gams_strict$p.fdr)## [1] 0.02839014nested_gams <- nested_gams %>%
mutate(LB = map(data, ~round(min(na.omit(.x$value)))),
UB = map(data, ~round(max(na.omit(.x$value)))),
range = map2(LB, UB, ~seq(from=.x, to=.y, by=1)),
logpredict = map2(model, range, ~predict(.x, newdata = data.frame(value = .y), type="link",se.fit=TRUE)),
fit = map(logpredict, ~plogis(.x$fit)),
lCI = map(logpredict, ~plogis(.x$fit-1.96*.x$se.fit)),
hCI = map(logpredict, ~plogis(.x$fit+1.96*.x$se.fit)))Define theme for gam plots
2a.1 Probability of exclusion as a function of ADOS
ados <- nested_gams %>%
filter(variable=="ADOS") %>%
select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>%
unnest(c(range, fit, lCI, hCI))
p_ados <- ggplot(ados, aes(x=range, y=fit))+
geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+
scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
labs(x='', y='Probability of Exclusion', fill='Motion Control', colour='Motion Control')+
scale_x_continuous(expand = c(0, 0))+
gam_theme+
ggtitle("ADOS")+
theme(plot.title = element_text(size = 11, hjust = 0.5))2a.2 Probability of exclusion as a function of SRS
srs <- nested_gams %>%
filter(variable=="SRS") %>%
select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>%
unnest(c(range, fit, lCI, hCI))
p_srs <- ggplot(srs, aes(x=range, y=fit))+
geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+
scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
labs(x='', y='', fill='Motion Control', colour='Motion Control')+
scale_x_continuous(expand = c(0, 0))+
gam_theme+
ggtitle("SRS")+
theme(plot.title = element_text(size = 11, hjust = 0.5))+
theme(axis.title.y = element_blank())2a.3 Probability of exclusion as a function of Inattention
ina <- nested_gams %>%
filter(variable=="Inattention") %>%
select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>%
unnest(c(range, fit, lCI, hCI))
p_in <- ggplot(ina, aes(x=range, y=fit))+
geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+
scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
labs(x='', y='', fill='Motion Control', colour='Motion Control')+
scale_x_continuous(expand = c(0, 0))+
gam_theme+
ggtitle("Inattention")+
theme(plot.title = element_text(size = 11, hjust = 0.5))+
theme(axis.title.y = element_blank())2a.4 Probability of exclusion as a function of Hyperactivity
hi <- nested_gams %>%
filter(variable=="Hyperactivity") %>%
select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>%
unnest(c(range, fit, lCI, hCI))
p_hi <- ggplot(hi, aes(x=range, y=fit))+
geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+
scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
labs(x='', y='', fill='Motion Control', colour='Motion Control')+
scale_x_continuous(expand = c(0, 0))+
gam_theme+
ggtitle("Hyperactivity")+
theme(plot.title = element_text(size = 11, hjust = 0.5))+
theme(axis.title.y = element_blank())2a.5 Probability of exclusion as a function of Motor Overflow
mo <- nested_gams %>%
filter(variable=="Motor Overflow") %>%
select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>%
unnest(c(range, fit, lCI, hCI))
p_mo <- ggplot(mo, aes(x=range, y=fit))+
geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+
scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
labs(x='', y='', fill='Motion Control', colour='Motion Control')+
scale_x_continuous(expand = c(0, 0))+
gam_theme+
ggtitle("Motor Overflow")+
theme(plot.title = element_text(size = 11, hjust = 0.5))+
theme(axis.title.y = element_blank())2a.6 Probability of exclusion as a function of Age
age<- nested_gams %>%
filter(variable=="Age") %>%
select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>%
unnest(c(range, fit, lCI, hCI))
p_age <- ggplot(age, aes(x=range, y=fit))+
geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+
scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
labs(x='', y='', fill='Motion Control', colour='Motion Control')+
scale_x_continuous(expand = c(0, 0))+
gam_theme+
ggtitle("Age")+
theme(plot.title = element_text(size = 11, hjust = 0.5))+
theme(axis.title.y = element_blank())2a.7 Probability of exclusion as a function of GAI
gai <- nested_gams %>%
filter(variable=="GAI") %>%
select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>%
unnest(c(range, fit, lCI, hCI))
p_gai <- ggplot(gai, aes(x=range, y=fit))+
geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+
scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
labs(x='', y='', fill='Motion Control', colour='Motion Control')+
scale_x_continuous(expand = c(0, 0))+
gam_theme+
ggtitle("GAI")+
theme(plot.title = element_text(size = 11, hjust = 0.5))+
theme(axis.title.y = element_blank())
p_legend = cowplot::get_legend(p_gai + guides(color = guide_legend(nrow = 1))+theme(legend.position = "bottom", legend.text = element_text(size = 10), legend.key.size=unit(.1, "in")))2b. Density plots of covariates used to fit GAMs (across included & excluded children)
Define theme for density plots of covariates across included and excluded children
2b.1 ADOS density
ddata <- nested_gams %>%
filter(variable=="ADOS") %>%
filter(Motion.Exclusion.Level=="Lenient") %>%
select("variable", "Motion.Exclusion.Level", "data") %>%
unnest(data) %>%
filter(PrimaryDiagnosis=="ASD")
ddata$PrimaryDiagnosis <- droplevels(ddata$PrimaryDiagnosis)
d_ados=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+
geom_density(alpha=0.5, inherit.aes=TRUE)+
scale_x_continuous(expand = c(0, 0))+
scale_y_continuous(expand = c(0, 0), limits = c(0, .09), breaks=seq(0, .08, by=.02))+
labs(x='', y='Density')+
scale_fill_manual(values = c("#FDE599"))+
scale_color_manual(values = c("#E9D38D"))+
den_theme2b.2 SRS density
ddata <- nested_gams %>%
filter(variable=="SRS") %>%
filter(Motion.Exclusion.Level=="Lenient") %>%
select("variable", "Motion.Exclusion.Level", "data") %>%
unnest(data)
d_srs=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+
geom_density(alpha=0.5, inherit.aes=TRUE)+
scale_x_continuous(expand = c(0, 0), limits=c(0,max(srs$range)),breaks = seq(0, 100 , by = 50))+
scale_y_continuous(expand = c(0, 0))+
scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), panel.background = element_blank(),
axis.line = element_line(size=.5), panel.border = element_blank())+
den_theme+
theme(axis.title.y = element_blank())+
labs(x='')2b.3 Inattention density
ddata <- nested_gams %>%
filter(variable=="Inattention") %>%
filter(Motion.Exclusion.Level=="Lenient") %>%
select("variable", "Motion.Exclusion.Level", "data") %>%
unnest(data)
d_in=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+
geom_density(alpha=0.5, inherit.aes=TRUE)+
scale_x_continuous(expand = c(0, 0))+
scale_y_continuous(expand = c(0, 0))+
scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
labs(x='')+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), panel.background = element_blank(),
axis.line = element_line(size=.5), panel.border = element_blank())+
theme(axis.title.y = element_blank())+
den_theme+
theme(axis.title.y = element_blank())+
labs(x='')2b.4 Hyperactivity/Impulsivity Density
ddata <- nested_gams %>%
filter(variable=="Hyperactivity") %>%
filter(Motion.Exclusion.Level=="Lenient") %>%
select("variable", "Motion.Exclusion.Level", "data") %>%
unnest(data)
d_hi=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+
geom_density(alpha=0.5, inherit.aes=TRUE)+
scale_x_continuous(expand = c(0, 0))+
scale_y_continuous(expand = c(0, 0))+
scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), panel.background = element_blank(),
axis.line = element_line(size=.5), panel.border = element_blank())+
den_theme+
theme(axis.title.y = element_blank())+
labs(x='')2b.5 Motor Overflow Density
ddata <- nested_gams %>%
filter(variable=="Motor Overflow") %>%
filter(Motion.Exclusion.Level=="Lenient") %>%
select("variable", "Motion.Exclusion.Level", "data") %>%
unnest(data)
d_mo=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+
geom_density(alpha=0.5, inherit.aes=TRUE)+
scale_x_continuous(expand = c(0, 0))+
scale_y_continuous(expand = c(0, 0))+
theme(axis.title.y = element_blank())+
scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), panel.background = element_blank(),
axis.line = element_line(size=.5), panel.border = element_blank())+
den_theme+
theme(axis.title.y = element_blank())+
labs(x='')2b.6 Age Density
ddata <- nested_gams %>%
filter(variable=="Age") %>%
filter(Motion.Exclusion.Level=="Lenient") %>%
select("variable", "Motion.Exclusion.Level", "data") %>%
unnest(data)
d_age=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+
geom_density(alpha=0.5, inherit.aes=TRUE)+
scale_x_continuous(expand = c(0, 0), limits=c(8,13), breaks = seq(8, 13 , by = 1))+
scale_y_continuous(expand = c(0, 0), limits=c(0,.29), breaks=seq(0, .25, by = .05))+
scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), panel.background = element_blank(),
axis.line = element_line(size=.5), panel.border = element_blank())+
den_theme+
theme(axis.title.y = element_blank())+
labs(x='')2b.7 GAI Density
ddata <- nested_gams %>%
filter(variable=="GAI") %>%
filter(Motion.Exclusion.Level=="Lenient") %>%
select("variable", "Motion.Exclusion.Level", "data") %>%
unnest(data)
d_gai=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+
geom_density(alpha=0.5, inherit.aes=TRUE)+
scale_x_continuous(expand = c(0, 0))+
scale_y_continuous(expand = c(0, 0), limits = c(0, .035), breaks=seq(0., .03, by=.01))+
scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), panel.background = element_blank(),
axis.line = element_line(size=.5), panel.border = element_blank())+
den_theme+
labs(x='')+
theme(axis.title.y = element_blank())
hist_legend = cowplot::get_legend(d_gai + guides(color = guide_legend(nrow = 1))+theme(legend.position = "bottom", legend.text = element_text(size = 10), legend.key.size=unit(.1, "in")))combine gam plots with densities & print
top_row <- cowplot::plot_grid(p_ados, p_srs, p_in, p_hi, p_mo, p_age, p_gai, ncol=7, rel_widths=c(1.18/7, .97/7, .97/7, .97/7, .97/7, .97/7, .97/7))
bottom_row <- cowplot::plot_grid(d_ados, d_srs, d_in, d_hi, d_mo, d_age, d_gai, ncol=7, rel_widths=c(1.18/7, .97/7, .97/7, .97/7, .97/7, .97/7, .97/7))## Warning: Removed 91 rows containing non-finite values (stat_density).## Warning: Removed 19 rows containing non-finite values (stat_density).png("Figures/fig_probExclusion_allGAM_TD_ASD_cc.png",width=10,height=6,units="in",res=200)
cowplot::plot_grid(p_legend, top_row, NULL, bottom_row, NULL, hist_legend, nrow=6, rel_heights=c(.1, 1, -.07, .5, -.07, .1))
dev.off()## quartz_off_screen
## 2#png("~/Dropbox/Apps/Overleaf/MotionSelectionBias_rsfMRI/Figures/fig_probExclusion_allGAM_TD_ASD_cc.png",width=10,height=6,units="in",res=200)
cowplot::plot_grid(p_legend, top_row, NULL, bottom_row, NULL, hist_legend, nrow=6, rel_heights=c(.1, 1, -.07, .5, -.07, .1))
#dev.off()NOTE: 19 rows = # of participants with missing Motor Overflow scores (imputed for deconfounded group difference)
3. Impact of motion QC on phenotypic representation within diagnostic groups
3a. Split violin plots
My_Theme = theme_light()+theme(
legend.title = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(size = 5),
axis.text.y = element_text(size = 8),
strip.text.x = element_text(size = 12, face = "bold", color="black"),
strip.text.y = element_text(size = 10, color="black"),
strip.background = element_rect(fill="white"),
plot.title = element_text(size = 9, hjust = 0.5))3a.1 ADOS (ASD only) split violin
NOTE: Missing ADOS scores for one participant evaluated at CARD
ados <- aim1 %>%
filter(PrimaryDiagnosis=="ASD" & variable=="ADOS") %>%
dplyr::select(-c("PrimaryDiagnosis", "variable"))
ados_violin <- ggplot(ados, aes(Motion.Exclusion.Level, value, fill = Included, color = Included)) +
geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE, adjust = 1.5) +
stat_summary(fun = "mean", position = position_dodge(width = 0.5),
color="black", geom="point", aes(y=value))+
stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5),
color="black", geom="errorbar", width=.2)+
scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
My_Theme+
theme(legend.text = element_text(size = 8))+
theme(legend.position = "none")+
ggtitle("ADOS")3a.2 SRS split violin
srs <- filter(aim1G, variable=="SRS")
aim1_srs <- ggplot(srs, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
facet_grid(PrimaryDiagnosis~., scales = "fixed")+
stat_summary(fun = "mean", position = position_dodge(width = 0.5),
color="black", geom="point", aes(y=value))+
stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5),
color="black", geom="errorbar", width=.2)+
# geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
# geom_pointrange(
# data = summary_pheno,
# aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
# shape = 20,
# position = position_dodge(width = 0.9))+
scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
My_Theme+
theme(strip.text.y =element_blank())+
theme(legend.position = "none")+
ggtitle("SRS")
v_legend = cowplot::get_legend(aim1_srs + guides(color = guide_legend(nrow = 2))+
theme(legend.position = "left",
legend.text = element_text(size = 10),
legend.key.size=unit(.1, "in")))## Warning: Removed 273 rows containing non-finite values (stat_ydensity).## Warning: Removed 273 rows containing non-finite values (stat_summary).
## Warning: Removed 273 rows containing non-finite values (stat_summary).NOTE: 273/3 = 91, # of participants missing SRS
with(dat3, table(PrimaryDiagnosis, is.na(SRS.Score)))##
## PrimaryDiagnosis FALSE TRUE
## TD 263 90
## ASD 150 13a.3 Inattention Split Violin
inat <- filter(aim1G, variable=="Inattention")
paper_inat <- ggplot(inat, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
facet_grid(PrimaryDiagnosis~.)+
stat_summary(fun = "mean", position = position_dodge(width = 0.5),
color="black", geom="point", aes(y=value))+
stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5),
color="black", geom="errorbar", width=.2)+
# geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
#geom_pointrange(
# data = summary_pheno,
# aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
# shape = 20,
# position = position_dodge(width = 0.9))+
scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
My_Theme+
theme(strip.text.y = element_blank())+
theme(legend.position = "none")+
ggtitle("Inattention")3a.4 Hyperactivity/Impulsivity Spilt Violin
hyp <- filter(aim1G, variable=="Hyperactivity")
paper_hyp <- ggplot(hyp, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
facet_grid(PrimaryDiagnosis~.)+
stat_summary(fun = "mean", position = position_dodge(width = 0.5),
color="black", geom="point", aes(y=value))+
stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5),
color="black", geom="errorbar", width=.2)+
# geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
#geom_pointrange(
# data = summary_pheno,
# aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
# shape = 20,
# position = position_dodge(width = 0.9))+
scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
My_Theme+
theme(strip.text.y = element_blank())+
theme(legend.position = "none")+
ggtitle("Hyperactivity")3a.5 Motor Overflow Split Violin
overflow <- filter(aim1G, variable=="Motor Overflow")
aim1_of <- ggplot(overflow, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
facet_grid(PrimaryDiagnosis~.)+
stat_summary(fun = "mean", position = position_dodge(width = 0.5),
color="black", geom="point", aes(y=value))+
stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5),
color="black", geom="errorbar", width=.2)+
# geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
#geom_pointrange(
# data = summary_pheno,
# aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
# shape = 20,
# position = position_dodge(width = 0.9))+
scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
My_Theme+
theme(strip.text.y =element_blank())+
theme(legend.position = "none")+
ggtitle("Motor Overflow")3a.6 Age Split Violin
age <- filter(aim1G, variable=="Age")
aim1_age <- ggplot(age, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
facet_grid(PrimaryDiagnosis~.)+
stat_summary(fun = "mean", position = position_dodge(width = 0.5),
color="black", geom="point", aes(y=value))+
stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5),
color="black", geom="errorbar", width=.2)+
# geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
#geom_pointrange(
# data = summary_pheno,
# aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
# shape = 20,
# position = position_dodge(width = 0.9))+
scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
My_Theme+
theme(strip.text.y =element_blank())+
theme(legend.position = "none")+
ggtitle("Age")3a.7 GAI Split Violin
gai <- filter(aim1G, variable=="GAI")
aim1_gai <- ggplot(gai, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
facet_grid(PrimaryDiagnosis~.)+
stat_summary(fun = "mean", position = position_dodge(width = 0.5),
color="black", geom="point", aes(y=value))+
stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5),
color="black", geom="errorbar", width=.2)+
# geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
#geom_pointrange(
# data = summary_pheno,
# aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
# shape = 20,
# position = position_dodge(width = 0.9))+
scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
My_Theme+
theme(legend.position = "none")+
ggtitle("GAI")Combine split violins into one figure & save
3b. Mann-Whitney U tests to compare included vs excluded participants using lenient motion QC (13 tests)
#run lenient tests first
aim1tib <- tibble(filter(aim1, Motion.Exclusion.Level=="Lenient"))
aim1MW <- aim1tib %>%
group_by(PrimaryDiagnosis, variable) %>%
tidyr::nest()
#hypothesis: included children will have less severe symptoms. NOTE: ADOS only collected in ASD (9 tests)
nested_mw_less <- aim1MW %>%
filter(variable %in% c("SRS", "Inattention", "Hyperactivity",
"Motor Overflow")|(variable=="ADOS"&PrimaryDiagnosis=="ASD")) %>%
mutate(mwm = map(data, ~wilcox.test(value~Included, alternative="less", data = na.omit(.x))),
includedMean = map(data, ~mean(.x$value, data=filter(.x, Included=="Included"))),
excludedMean = map(data, ~mean(.x$value, data=filter(.x, Included=="Excluded"))),
coefs = map(mwm, tidy)) %>%
unnest(c(coefs, includedMean, excludedMean))
#hypothesis: included children will be older and have greater GAI (4 tests)
nested_mw_greater<- aim1MW %>%
filter(variable %in% c("Age", "GAI")) %>%
mutate(mwm = map(data, ~wilcox.test(value~Included, alternative="greater", data = na.omit(.x))),
includedMean = map(data, ~mean(.x$value, data=filter(.x, Included=="Included"))),
excludedMean = map(data, ~mean(.x$value, data=filter(.x, Included=="Excluded"))),
coefs = map(mwm, tidy)) %>%
unnest(c(coefs, includedMean, excludedMean))
complete_mw <- rbind(nested_mw_less, nested_mw_greater)
complete_mw$p.fdr <- p.adjust(complete_mw$p.value, method = "BH")
complete_mw[order(complete_mw$PrimaryDiagnosis, decreasing=TRUE), c(1:2, 5:9)]## # A tibble: 13 × 7
## # Groups: PrimaryDiagnosis, variable [13]
## PrimaryDiagnosis variable includedMean excludedMean statistic p.value method
## <fct> <fct> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 ASD ADOS 14.5 14.5 1780. 9.21e-3 Wilco…
## 2 ASD SRS NA NA 1720 5.79e-3 Wilco…
## 3 ASD Motor Ov… NA NA 1260 9.48e-4 Wilco…
## 4 ASD Inattent… 16.7 16.7 2442. 6.42e-1 Wilco…
## 5 ASD Hyperact… 12.0 12.0 2419 6.06e-1 Wilco…
## 6 ASD Age 10.3 10.3 2848 2.16e-2 Wilco…
## 7 ASD GAI 107. 107. 2960 6.55e-3 Wilco…
## 8 TD SRS NA NA 4562. 4.31e-1 Wilco…
## 9 TD Motor Ov… NA NA 7086. 5.69e-2 Wilco…
## 10 TD Inattent… 3.10 3.10 8004 2.68e-1 Wilco…
## 11 TD Hyperact… 2.13 2.13 7018. 1.98e-2 Wilco…
## 12 TD Age 10.4 10.4 10074. 1.02e-2 Wilco…
## 13 TD GAI 116. 116. 9882 2.02e-2 Wilco…3c. Mann-Whitney U tests to compare included vs excluded participants using strict motion QC (13 tests)
#run tests using strict motion QC
aim1tib <- tibble(filter(aim1, Motion.Exclusion.Level=="Strict"))
aim1MW <- aim1tib %>%
group_by(variable, PrimaryDiagnosis) %>%
tidyr::nest()
#hypothesis: included children will have less severe symptoms
nested_mw_less <- aim1MW %>%
filter(variable %in% c("SRS", "Inattention", "Hyperactivity",
"Motor Overflow")|(variable=="ADOS"&PrimaryDiagnosis=="ASD")) %>%
mutate(mwm = map(data, ~wilcox.test(value~Included, alternative="less", data = na.omit(.x))),
includedMean = map(data, ~mean(.x$value, data=filter(.x, Included=="Included"))),
excludedMean = map(data, ~mean(.x$value, data=filter(.x, Included=="Excluded"))),
coefs = map(mwm, tidy)) %>%
unnest(c(coefs, includedMean, excludedMean))
#hypothesis: included children will be older and have greater GAI
nested_mw_greater<- aim1MW %>%
filter(variable %in% c("Age", "GAI")) %>%
mutate(mwm = map(data, ~wilcox.test(value~Included, alternative="greater", data = na.omit(.x))),
includedMean = map(data, ~mean(.x$value, data=filter(.x, Included=="Included"))),
excludedMean = map(data, ~mean(.x$value, data=filter(.x, Included=="Excluded"))),
coefs = map(mwm, tidy)) %>%
unnest(c(coefs, includedMean, excludedMean))
complete_mw <- rbind(nested_mw_less, nested_mw_greater)
complete_mw$p.fdr <- p.adjust(complete_mw$p.value, method = "BH")
complete_mw[order(complete_mw$PrimaryDiagnosis, decreasing=TRUE), c(1:2, 5:6, 10:11)]## # A tibble: 13 × 6
## # Groups: PrimaryDiagnosis, variable [13]
## PrimaryDiagnosis variable includedMean excludedMean alternative p.fdr
## <fct> <fct> <dbl> <dbl> <chr> <dbl>
## 1 ASD ADOS 14.5 14.5 less 0.0767
## 2 ASD SRS NA NA less 0.0763
## 3 ASD Motor Overflow NA NA less 0.0847
## 4 ASD Inattention 16.7 16.7 less 0.223
## 5 ASD Hyperactivity 12.0 12.0 less 0.322
## 6 ASD Age 10.3 10.3 greater 0.0763
## 7 ASD GAI 107. 107. greater 0.303
## 8 TD SRS NA NA less 0.0847
## 9 TD Motor Overflow NA NA less 0.194
## 10 TD Inattention 3.10 3.10 less 0.194
## 11 TD Hyperactivity 2.13 2.13 less 0.0159
## 12 TD Age 10.4 10.4 greater 0.0847
## 13 TD GAI 116. 116. greater 0.111xtable(complete_mw[order(complete_mw$PrimaryDiagnosis, decreasing=TRUE), c(1:2, 5:6, 10:11)])## % latex table generated in R 4.1.1 by xtable 1.8-4 package
## % Thu Oct 7 00:44:19 2021
## \begin{table}[ht]
## \centering
## \begin{tabular}{rllrrlr}
## \hline
## & PrimaryDiagnosis & variable & includedMean & excludedMean & alternative & p.fdr \\
## \hline
## 1 & ASD & ADOS & 14.54 & 14.54 & less & 0.08 \\
## 2 & ASD & SRS & & & less & 0.08 \\
## 3 & ASD & Motor Overflow & & & less & 0.08 \\
## 4 & ASD & Inattention & 16.67 & 16.67 & less & 0.22 \\
## 5 & ASD & Hyperactivity & 11.99 & 11.99 & less & 0.32 \\
## 6 & ASD & Age & 10.32 & 10.32 & greater & 0.08 \\
## 7 & ASD & GAI & 107.11 & 107.11 & greater & 0.30 \\
## 8 & TD & SRS & & & less & 0.08 \\
## 9 & TD & Motor Overflow & & & less & 0.19 \\
## 10 & TD & Inattention & 3.10 & 3.10 & less & 0.19 \\
## 11 & TD & Hyperactivity & 2.13 & 2.13 & less & 0.02 \\
## 12 & TD & Age & 10.36 & 10.36 & greater & 0.08 \\
## 13 & TD & GAI & 115.54 & 115.54 & greater & 0.11 \\
## \hline
## \end{tabular}
## \end{table}4. Edgewise functional connectivity as a function of the covariates
Combine partial correlations with melted rs-fMRI usability and covariate info
startEdgeidx = which(names(dat3)=='r.ic1.ic2')
endEdgeidx = which(names(dat3)=='r.ic29.ic30')
names(dat3)[startEdgeidx:endEdgeidx]## [1] "r.ic1.ic2" "r.ic1.ic4" "r.ic1.ic8" "r.ic1.ic13" "r.ic1.ic14"
## [6] "r.ic1.ic15" "r.ic1.ic17" "r.ic1.ic19" "r.ic1.ic21" "r.ic1.ic22"
## [11] "r.ic1.ic24" "r.ic1.ic25" "r.ic1.ic26" "r.ic1.ic27" "r.ic1.ic28"
## [16] "r.ic1.ic29" "r.ic1.ic30" "r.ic2.ic4" "r.ic2.ic8" "r.ic2.ic13"
## [21] "r.ic2.ic14" "r.ic2.ic15" "r.ic2.ic17" "r.ic2.ic19" "r.ic2.ic21"
## [26] "r.ic2.ic22" "r.ic2.ic24" "r.ic2.ic25" "r.ic2.ic26" "r.ic2.ic27"
## [31] "r.ic2.ic28" "r.ic2.ic29" "r.ic2.ic30" "r.ic4.ic8" "r.ic4.ic13"
## [36] "r.ic4.ic14" "r.ic4.ic15" "r.ic4.ic17" "r.ic4.ic19" "r.ic4.ic21"
## [41] "r.ic4.ic22" "r.ic4.ic24" "r.ic4.ic25" "r.ic4.ic26" "r.ic4.ic27"
## [46] "r.ic4.ic28" "r.ic4.ic29" "r.ic4.ic30" "r.ic8.ic13" "r.ic8.ic14"
## [51] "r.ic8.ic15" "r.ic8.ic17" "r.ic8.ic19" "r.ic8.ic21" "r.ic8.ic22"
## [56] "r.ic8.ic24" "r.ic8.ic25" "r.ic8.ic26" "r.ic8.ic27" "r.ic8.ic28"
## [61] "r.ic8.ic29" "r.ic8.ic30" "r.ic13.ic14" "r.ic13.ic15" "r.ic13.ic17"
## [66] "r.ic13.ic19" "r.ic13.ic21" "r.ic13.ic22" "r.ic13.ic24" "r.ic13.ic25"
## [71] "r.ic13.ic26" "r.ic13.ic27" "r.ic13.ic28" "r.ic13.ic29" "r.ic13.ic30"
## [76] "r.ic14.ic15" "r.ic14.ic17" "r.ic14.ic19" "r.ic14.ic21" "r.ic14.ic22"
## [81] "r.ic14.ic24" "r.ic14.ic25" "r.ic14.ic26" "r.ic14.ic27" "r.ic14.ic28"
## [86] "r.ic14.ic29" "r.ic14.ic30" "r.ic15.ic17" "r.ic15.ic19" "r.ic15.ic21"
## [91] "r.ic15.ic22" "r.ic15.ic24" "r.ic15.ic25" "r.ic15.ic26" "r.ic15.ic27"
## [96] "r.ic15.ic28" "r.ic15.ic29" "r.ic15.ic30" "r.ic17.ic19" "r.ic17.ic21"
## [101] "r.ic17.ic22" "r.ic17.ic24" "r.ic17.ic25" "r.ic17.ic26" "r.ic17.ic27"
## [106] "r.ic17.ic28" "r.ic17.ic29" "r.ic17.ic30" "r.ic19.ic21" "r.ic19.ic22"
## [111] "r.ic19.ic24" "r.ic19.ic25" "r.ic19.ic26" "r.ic19.ic27" "r.ic19.ic28"
## [116] "r.ic19.ic29" "r.ic19.ic30" "r.ic21.ic22" "r.ic21.ic24" "r.ic21.ic25"
## [121] "r.ic21.ic26" "r.ic21.ic27" "r.ic21.ic28" "r.ic21.ic29" "r.ic21.ic30"
## [126] "r.ic22.ic24" "r.ic22.ic25" "r.ic22.ic26" "r.ic22.ic27" "r.ic22.ic28"
## [131] "r.ic22.ic29" "r.ic22.ic30" "r.ic24.ic25" "r.ic24.ic26" "r.ic24.ic27"
## [136] "r.ic24.ic28" "r.ic24.ic29" "r.ic24.ic30" "r.ic25.ic26" "r.ic25.ic27"
## [141] "r.ic25.ic28" "r.ic25.ic29" "r.ic25.ic30" "r.ic26.ic27" "r.ic26.ic28"
## [146] "r.ic26.ic29" "r.ic26.ic30" "r.ic27.ic28" "r.ic27.ic29" "r.ic27.ic30"
## [151] "r.ic28.ic29" "r.ic28.ic30" "r.ic29.ic30"signalFC <- dat3[, c(1,startEdgeidx:endEdgeidx)]
dat2 <- merge(aim1, signalFC, all=TRUE, by = "ID")
fcMelt <- reshape2::melt(dat2[,1:160],
id.vars=names(dat2)[1:7],
variable.name = "edge",
value.name = "fc")4a. Run nested gams
fctib <- tibble(filter(fcMelt, Motion.Exclusion.Level!="None"))
fctib$Motion.Exclusion.Level <- droplevels(fctib$Motion.Exclusion.Level)
fcNest <- fctib %>%
filter(Included=="Included") %>%
group_by(variable, Motion.Exclusion.Level, edge) %>%
tidyr::nest()
#nested models
fc_gams <- fcNest %>%
mutate(model = map(data, ~mgcv::gam(fc~s(value), data = na.omit(.x), method="REML")),
coefs = map(model, tidy, conf.int = FALSE)) %>%
unnest(coefs)4b. Plot histograms of edgewise p-values from GAMs
4b.1 Histogram of edgewise p-values for partial correlations as a function of ADOS.
NOTE: TD scores = 0
pdat <- fc_gams %>%
ungroup() %>%
filter(variable=="ADOS") %>%
select(Motion.Exclusion.Level, p.value) %>%
group_by(Motion.Exclusion.Level)
hist_ados_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+
geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
scale_x_continuous(expand = c(0, 0))+
scale_y_continuous(expand = c(0, 0))+
scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), panel.background = element_blank(),
axis.line = element_line(size=.5), panel.border = element_blank())+
My_Theme+
ggtitle("ADOS")+
theme(plot.title = element_text(size = 9, hjust = 0.5))+
theme(legend.position = "none")+
labs(x='p value', y='Count')+
theme(axis.title.y = element_text(size = 9, angle=90))+
theme(axis.title.x = element_text(size = 7))4b.2 Histogram of edgewise p-values for partial correlations as a function of SRS
pdat <- fc_gams %>%
ungroup() %>%
filter(variable=="SRS") %>%
select(Motion.Exclusion.Level, p.value) %>%
group_by(Motion.Exclusion.Level)
hist_srs_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+
geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
scale_x_continuous(expand = c(0, 0))+
scale_y_continuous(expand = c(0, 0))+
labs(x='p value', y='')+
scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), panel.background = element_blank(),
axis.line = element_line(size=.5), panel.border = element_blank())+
My_Theme+
theme(axis.title.y = element_blank())+
ggtitle("SRS")+
theme(plot.title = element_text(size = 9, hjust = 0.5))+
theme(legend.position = "none")+
theme(axis.title.x = element_text(size = 7))4b.3 Histogram of p-values for partial correlations as a function of inattention
pdat <- fc_gams %>%
ungroup() %>%
filter(variable=="Inattention") %>%
select(Motion.Exclusion.Level, p.value) %>%
group_by(Motion.Exclusion.Level)
hist_in_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+
geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
scale_x_continuous(expand = c(0, 0))+
scale_y_continuous(expand = c(0, 0))+
labs(x='p value', y='')+
scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), panel.background = element_blank(),
axis.line = element_line(size=.5), panel.border = element_blank())+
My_Theme+
theme(axis.title.y = element_blank())+
ggtitle("Inattention")+
theme(plot.title = element_text(size = 9, hjust = 0.5))+
theme(legend.position = "none")+
theme(axis.title.x = element_text(size = 7))4b.4 Histogram of p-values for partial correlations as a function of hyperactivity
pdat <- fc_gams %>%
ungroup() %>%
filter(variable=="Hyperactivity") %>%
select(Motion.Exclusion.Level, p.value) %>%
group_by(Motion.Exclusion.Level)
hist_hi_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+
geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
scale_x_continuous(expand = c(0, 0))+
scale_y_continuous(expand = c(0, 0))+
labs(x='p value', y='')+
scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), panel.background = element_blank(),
axis.line = element_line(size=.5), panel.border = element_blank())+
My_Theme+
theme(axis.title.y = element_blank())+
ggtitle("Hyperactivity")+
theme(plot.title = element_text(size = 9, hjust = 0.5))+
theme(legend.position = "none")+
theme(axis.title.x = element_text(size = 7))4b.5 Histogram of p-values for partial correlations as a function of motor overflow
pdat <- fc_gams %>%
ungroup() %>%
filter(variable=="Motor Overflow") %>%
select(Motion.Exclusion.Level, p.value) %>%
group_by(Motion.Exclusion.Level)
hist_mo_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+
geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
scale_x_continuous(expand = c(0, 0))+
scale_y_continuous(expand = c(0, 0))+
labs(x='p value', y='')+
scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), panel.background = element_blank(),
axis.line = element_line(size=.5), panel.border = element_blank())+
My_Theme+
theme(axis.title.y = element_blank())+
ggtitle("Motor Overflow")+
theme(plot.title = element_text(size = 9, hjust = 0.5))+
theme(legend.position = "none")+
theme(axis.title.x = element_text(size = 7))4b.6 Histogram of p-values for partial correlations as a function of age
pdat <- fc_gams %>%
ungroup() %>%
filter(variable=="Age") %>%
select(Motion.Exclusion.Level, p.value) %>%
group_by(Motion.Exclusion.Level)
hist_age_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+
geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
scale_x_continuous(expand = c(0, 0))+
scale_y_continuous(expand = c(0, 0))+
labs(x='p value', y='')+
scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), panel.background = element_blank(),
axis.line = element_line(size=.5), panel.border = element_blank())+
My_Theme+
theme(axis.title.y = element_blank())+
ggtitle("Age")+
theme(plot.title = element_text(size = 9, hjust = 0.5))+
theme(legend.position = "none")+
theme(axis.title.x = element_text(size = 7))4b.7 Histogram of p-values for partial correlations as a function of GAI
pdat <- fc_gams %>%
ungroup() %>%
filter(variable=="GAI") %>%
select(Motion.Exclusion.Level, p.value) %>%
group_by(Motion.Exclusion.Level)
hist_gai_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+
geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
scale_x_continuous(expand = c(0, 0))+
scale_y_continuous(expand = c(0, 0))+
labs(x='p value', y='')+
scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), panel.background = element_blank(),
axis.line = element_line(size=.5), panel.border = element_blank())+
My_Theme+
theme(axis.title.y = element_blank())+
ggtitle("GAI")+
theme(plot.title = element_text(size = 9, hjust = 0.5))+
theme(legend.position = "none")+
theme(axis.title.x = element_text(size = 7))
hist_p_legend <- cowplot::get_legend(hist_gai_p + guides(color = guide_legend(nrow = 1))+theme(legend.position = "bottom", legend.text = element_text(size = 10), legend.key.size=unit(.1, "in"), legend.title = element_blank()))4b.8 Combine histograms and print
fc_hist <- cowplot::plot_grid(hist_ados_p, hist_srs_p, hist_in_p, hist_hi_p, hist_mo_p, hist_age_p, hist_gai_p, ncol=7, rel_widths=c(1.18/7, .97/7, .97/7, .97/7, .97/7, .97/7, .97/7))
png("Figures/fig_hist_rfc_cc.png",width=8,height=3,units="in",res=200)
cowplot::plot_grid(fc_hist, hist_p_legend, nrow=2, rel_heights=c(1, .1))
dev.off()## quartz_off_screen
## 2cowplot::plot_grid(fc_hist, hist_p_legend, nrow=2, rel_heights=c(1, .1))